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On the basis of the quasiclassical theory of superconductivity, we obtain a formula for the 
local density of states (LDOS) around a vortex core of superconductors with anisotropic pair- 
potential and Fermi surface in arbitrary directions of magnetic fields. Earlier results on the 
LDOS of d-wave superconductors and NbSe2 are naturally interpreted within our theory geo- 
metrically; the region with high intensity of the LDOS observed in numerical calculations turns 
out to the enveloping curve of the trajectory of Andreev bound states. We discuss experimental 
results on YNi2B2C within the quasiclassical theory of superconductivity. 
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1. Introduction 

One of the aims of the study on the condensed matter 
is to identify the order parameter in a symmetry-broken 
state. In supercoducting materials, the order parameter 
is given by the pairing wavefunction, which is related 
to the pair-potential. As experimental techniques to ob- 
serve the properties of the pair-potential, several probes 
are available such as nuclear-spin-lattice relaxation, the 
measurement of specific heat, (angular-resolved) photo- 
emission spectroscopy, field direction dependence of ther- 
mal conductivity and so on. Among those probes, we 
focus particularly on the scanning tunneling microscopy 
and spectroscopy (STM/STS), which is now a main tech- 
nology in nano-science. 

The observation of the local density of states (LDOS) 
in the vortex phase in superconductors yields rich in- 
formation on the pairing symmetry; the LDOS around 
a vortex core mainly comes from Andreev bound states, 
whose spectrum is sensitive to the anisotropy of the pair- 
potential. Thus the anisotropic spatial pattern of the 
LDOS around a vortex core is expected to be observed 
as a consequence of the anisotropy of pair-potential on 
the Fermi surface. Indeed the LDOS of NbSe2 1 and 
YM2B2C 2 measured by STM reveal, respectively, six- 
fold star-shaped and fourfold structure. Motivated by the 
experimental results on NbSe2, Gygi and Schliiter 3 and 
Hayashi et al. A successfully reproduced the LDOS around 
a vortex by numerical calculations. The LDOS around 
a vortex in non-s-wave superconductors was firstly ob- 
tained by Schopohl and Maki 5 in a numerical study of 
(i-wave vortex with the use of Riccati formalism of qua- 
siclassical theory. Ichioka et a/., 6 subsequently, obtained 
the LDOS of d-wave vortex in the self-consistent numer- 
ical calculations. Prior to ref. 5, Riccati formalism was 
developed by Ashida et al. 7 in the context of bound- 
ary problem at Supercondutor- Normal- metal interfaces. 8 
The formalisms of refs. 7 and 9 turned out to be equiva- 



lent although it is not so obvious. Riccati formalism was 
extended to non-equilibrium case by Eschrig 10 via the 
projection operator formalism by Shelankov. 11 

Although there exist many numerical works on the 
LDOS, each calculation usually takes much machine- 
time. It would be elaborate to treat three-dimensional 
superconductors with anisotropic Fermi surface and pair- 
ing interactions. Further, the LDOS of refs. 4 and 5 has 
strong spatial dependence and hence a comprehensive in- 
terpretation by an analytical theory is highly desirable. 
The bound states localized near an isotropic s-wave vor- 
tex core were studied analytically by Caroli et al 12 with 
the use of Bogoliubov-de Gennes equation. Equivalently, 
an approximate expression for the quasiclassical Green 
function near a vortex core for low energy was analyt- 
ically obtained by Kramer and Pesch. 13 We call their 
approximation "Kramer-Pesch approximation (KPA)" in 
the present paper. Although the spirit of KPA is hard to 
understand as it stands, it becomes clear after transform- 
ing quasiclassical equation into Riccati formalism; KPA 
is equivalent to first-order perturbation with respect to 
energy and impact parameter. KPA was successfully ap- 
plied to two-dimensional <i-wave superconductors, 5 in a 
slightly different formalism. 14 

In this paper, we apply KPA to more general case; 
the LDOS of two-dimensional and three-dimensional su- 
perconductors with anisotropic Fermi surface and pair- 
potential for arbitrary direction of magnetic fields are 
discussed. Our theory yields a natural geometrical inter- 
pretation of spatially anisotropic pattern of the LDOS for 
anisotropic superconductors. We will see that earlier nu- 
merical results on the LDOS are easily understood within 
our theory. Furthermore, we discuss experimental results 
on YM2B2C. 2 Contrary to the claim of the authors of 
ref. 2, we show that a part of the experimental results 
is naturally interpreted by the quasiclassical theory. In 
an earlier short report, the application of our theory to 
CePtsSi was already reported. 15 One of the aims of this 
paper is to present the details of the derivation. The un- 
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published results of ref. 16 on two-dimensional supercon- 
ductors are contained as a part of this paper (§3.1.2). 

This paper is organized as follows. In §2, we summarize 
the formulation of the quasiclassical theory of supercon- 
ductivity. In §3, we present our analytical formulation for 
the LDOS around a vortex core in the case of the unitary 
pair-potential (AA* oc &q) of the spin singlet or the spin 
triplet, using the Riccati formalism and Kramer-Pesch 
approximation. In §4, we calculate the specific LDOS 
patterns about NbSe2 and YM2B2C with our formula- 
tion. In §5, some remarks on the formulation are given. 
In §6, we compare our results with those of STM mea- 
surements on YM2B2C. In §7, the summary is given. In 
the Appendix A, we describe the formulation to obtain 
the matrix Riccati equation. In the Appendix B, we cal- 
culate the LDOS of CePtsSi as an example of the case 
of the pair-potential which does not satisfy the relation 
A A 1 " oc <r . 

2. Quasiclassical Theory of Superconductivity 

We consider the pair-potential: 

A m ,(R,r) = {^(fl+rjVy (i*-0),(2.1) 



A m ,(R, k) 



I 



dre- ikr 'V(ii,r), 



(2.2) 



where r and R are the relative and center-of-mass cood- 
inates, respectively. Here, ip v ( r ) is the field operator for 
the quasiparticle state with spin n and we use units in 
which h = &b = 1- We assume the weak coupling in- 
teraction so that the pair-potential is not zero only near 
the Fermi surface. We consider pair-potential written as 
either A(R, k) = iip(k)a y A(R) for singlet pairing or 
A(R, k) = i(d(k) • &)a y A(R) for triplet pairing. Here 
A(R) is a function of R. If the coherence length £ is 
large compared to the Fermi wave length (~ 1/&f), we 
can calculate the LDOS around a vortex core on the ba- 
sis of the quasiclassical theory of superconductivity. 17 19 
We consider the quasiclassical Green function g that has 
the matrix elements in the Nambu (particle-hole) space 



matrix. Considering clean superconductors in the type 
II limit, we neglect the impurity effect and the vector 
potential. 

The local density of states with the isotropic Fermi 
surface is given by 

Kr,6) = -K0) / ^Retr(5 R ). (2.6) 

Here, g R is retarded Green function and v(0) denotes the 
density of states on Fermi surface in the normal metalic 
state. The local density of states with the anisotropic 
Fermi surface is given by 



Here, dSp is the Fermi-surface area element and v F is the 
modulus of Fermi velocity. 

3. Riccati Formalism and Kramer-Pesch Ap- 
proximation 

In order to solve the quasiclassical equation (2.4), we 
simplify the equation with a parametrization for the 
propagators that satisfy the normalization condition. 
Projection operators are defined by P± = ^(1 =p g). 10 ' n 
Using this projection method (See, Appendix A), we ob- 
tain matrix Riccati equations: 



vf • Va + 2uo n a + aA^a — A 
v F • Vb - 2uj n b - bAb + A 1 " 



0, 
0, 



where 



-N 



(1 - ab) 
-2ib 



N 



(1+aS)" 1 




2ia 
-(i - ba) 



(1 + Sa)" 1 



(3.8) 
(3.9) 

(3.10) 
(3.11) 



Since these matrix Riccati equations (3.8) and (3.9) con- 
tain V only through v F • V, they reduce to a one- 
dimensional problem on a straight line, the direction of 
which is given by that of the Fermi velocity. 



g(r,k,\u n ) 




(2.3) 



where uo n = 7rT(2n + 1) is the Matsubara frequency with 
the temperature T and an integer n. Here k is a unit 
vector (fc = k/kp). Throughout the paper, "hat" a de- 
notes the 2x2 matrix in the spin space, and "check" 
a denotes a 4 x 4 matrix composed of the 2x2 Nambu 
space and the 2x2 spin space. 20 The equation of motion 
for g called Eilenberger equation is written as 17 ' 18 



iv F (fe) • V<? + [kj n f 3 - A, g] = 0, 



where 



T3 = 



<to 
-ao 



A = 



A 
-At 



(2.4) 



(2.5) 



Here vp(k) is the Fermi velocity. The commutator 
[a, b] = ab — ba. The Green function g satisfies the nor- 
malization condition: g 2 = 1, where 1 is the 4x4 unit 



3.1 Two-dimensional Fermi surface 
3.1.1 Frame of the space 

We consider the system with two-dimensional Fermi 
surface. We introduce the frame of the space written as 



k 



v 
u 



cos 6a + 4 / — - sin Ob 

m 




cos 6a + * / — sin 6b 

m 



v F (6 v ) 



cos 6 v a + sin 6 v b 



Xd + Yb, 
xv + yu, 



cos 6 V sin 6 V 
— sin 6 V cos 6 V 



(3.12) 
(3.13) 

(3.14) 

(3.15) 
(3.16) 

(3.17) 



Here, 6 V is the angle between the a axis and the Fermi ve- 
locity (See, Fig. 1) and m x and m y are the anisotropic ef- 
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fective masses. k a and fc& are the axes fixed to the crystal 
axes in momentum space. In this system, a point (k x ,k y ) 
on the Fermi surface with the energy e = k 2 /(2m x ) + 
k 2 /(2m y ) being the Fermi energy is parameterized with 
k x = kcosO and k y = k^/m y /m x sin#. 



4 

± b 




\ Y 


r 












\X a' 



Fig. 1. Frame of the space. The direction of v is equal to that of 
Fermi velocity, a and b are the axes fixed to the crystal axes. r 
is the angle between the a axis and r. 



3.1.2 Isotropic Fermi surface 21 

First, we assume superconductors with the circu- 
lar Fermi surface (m x = m y = m). Those with the 
anisotropic Fermi surface will be discussed briefly in 
§3.1.3. 

In the case of A At oc <to, we can solve the matrix 
Riccati equations easily. Equations (3.8) and (3.9) are to 
be solved, under the initial conditions: 



uo 2 n + ^trAAt 



At 



for x 



for i 



(3.18) 
(3.19) 



where Xq is a positive real cut-off. 

At the vortex core and low energy, the Green func- 
tion diverges, 13 which owes to the Andreev bound state. 
Therefore, we expand the denominator of the Green func- 
tion around a vortex core. We call this Kramer-Pesch ap- 
proximation, which is appropriate around a vortex core 
in a low energy region(|cc; n | <C |Aoo|). Here, Aoo is a 
pair-potential in the bulk region. We consider the per- 
turbation with respect to an impact parameter y and an 
energy uo n . We consider a pair-potential having the form 
of 



A(r,fc) 



A f (r,fe) 



f{r) exp(i# r )A(fe), 



(3.20) 



/(r)expK)^±i^A(/c), (3.21) 



V x<2 + y* 

/(r)exp(-i(9 r )At(fe), 



f(r) exp(-i0 v ) 



x-iy 



\A 2 +y 2 



(3.22) 
At (fe). (3.23) 



Here, f(r) — f(\/x 2 + y 2 ) is the spatial variation of pair- 
potential and /(0) = 0, lim r ^ 00 f(r) = A^ and r is an 
angle from the a axis (See, Fig. 1). We expand a pair- 
potential written as 

A(r, fe) ~ f(\x\) exp(i^) sign(x) (l + i|) A(fe),(3.24) 

= Ao + Ai, (3.25) 

At(r, fe) ~ f(\x\) exp(-i^) sign(x) (l - i|) At (fe), (3.26) 

= At+Al, (3.27) 

with respect to y. Here, Ai is defined as i(y/x)Ao. We 
should notice that Ai does not have the singularity at 
x = since f(\x\)/xis finite everywhere. Further remarks 
will be presented in §5. 

In the zero-th order with respect to y and uo ni equa- 
tions (3.8) and (3.9) are written as 



0, 



dbo 

dx 



L'F— - 6 A 6 + AJ = 0. 



(3.28) 
(3.29) 



Under the initial conditions (3.18) and (3.19), we can 
obtain the solutions: 



a 



bo 



A 



V^ trA o A o 



-A(fc)exp(i0 v ) 
AW ' 

At(fc)exp(-i0 v ) 



(3.30) 
(3.31) 



where A(fe) = y ^trA(fe)At(fe). For the moment, we de- 
note A(fe) by A for simplicity. 

In the first order with respect to y and uo ni equations 
(3.8) and (3.9) are written as 



° F ~dx + 2cJn ^° + ^oAJai + ^lAjao 



_|_ A . A t ; 

+a AJa - Ai 



0,(3.32) 



db 

v F — 2co n b - boAoh - 6iA 6 

ox 



-6 Ai6 + At = 0.(3.33) 

These equations are inhomogeneous linear differential 
equations. Therefore, the solutions 

2 



ai 



v F 



exp [2XF(x)] 



[ dx' 

J — x c 



y ~ 

i— Aq — cj n ao exp [— 2XF(x f 



+ exp [2A{ F(x) - F(-x c ) }] exp(i0 v ) 



(-1) 
A 



. y 



A^A 



A(fe), 



(3.34) 



h = 



v F 



exp [2AF(z)] 
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dx' 



y"0 



- uj n bo 



exp [—2XF(x') 



+ exp [2A{ F(x) - F(x c ) }] exp(-i^) 



A V x c . 



At(fe) 



(3.35) 



are obtained by the method of variation of constants. 
Here we introduce 



F(x) 



i 

- / dx'f(x f ). (3. 

V F JO 



36) 



Substituting a = do + ai and 6 = 6o + &i into eq. (3.10), 
we obtain Green function up to the first order: 



-2(a &i + aibo) x , 

VFe -2AF(.)^ o 



i^A/(i^i)+^ n 



e -2AF(V) ' 



-2XF(x) 



^0, 



2C [lE + c^ n ] 
in the limit x c — > oo. Here 



/>oo 

C(0) = / cte'e 
Jo 



-2AF(V) 



£7(5) 



c 



— 2AF(k') 



(3.37) 
,(3.38) 

(3.39) 

(3.40) 
(3.41) 



and E is independent of x. 

The retarded Green function g R is obtained by replac- 
ing \uo n with e + iS in eq. (3.39). When the Fermi surface 
is isotropic, V = (See eqs. (3.13) and (3.14)). Near the 
vortex core (|r| <C £), we use the following approxima- 
tions: 



v{r,e) 



lim 



1/(0) 



JO 



v F dO 5e~ 2XF W 



Co 



(3.42) 
Jo 



£^6 2 J 2ttC (e - £) 2 + S 2 
v(0)v F r ' 2 



where 



C 



E ■ 



2A Aqq^ 

e v f 



Jo 

yX /*° A^A 
Co Jo 



A' 



(3.44) 



dx 7 



y 



AooA 2 , (3.45) 



Co Jo Co Co 

with £o = v f /(ttA 00 ). Here, we note that the integrand 
of z/(r, e) in eq. (3.43) is not zero on the path where the 
relation (e — E(r, 0) = 0) is satisfied. The number of the 
solutions 0i satisfying e — E(r, 0i) = is finite as shown 
in Fig. 2. We call this path the quasiparticle path. 16 

We consider the points at which the LDOS diverges. 
We rewrite eq. (3.43) as 

U{0)VF d6 T A(6>) 6{6 ~ 6i) e -2AWF(x(0» 
L A^i y ' dh(r,e,e) ' 
* m 

(3.46) 

dh(r M exp{-2A(^) J F(a ; (^))} ) 

(3.47) 



f(r,e) 



Co 



z/(0)^f 



do 



with h(r,0,e) = e — E(r,0). The condition that the 




Fig. 2. Paths that contribute to the integrand of the LDOS at a 
certain point r in real space for an isotropic s-wave supercon- 
ductor. O denotes the vortex center. 6\ and 62 are the angles 
characterizing these paths. 



LDOS diverges at a given point r for an energy e is given 

by 



h(r,0,e) 

dh(r,0,e) 
d~0 



0, 



= 0. 



(3.48) 
(3.49) 



When we regard eq. (3.48) as the equations of a family of 
lines on x-y plane with one-parameter equations (3.48) 
and (3.49) describe the enveloping curve of that family 
of lines (= the quasiparticle paths in the present case). 
From eqs. (3.48) and (3.49), the parametric equation 



<9A(6>) 

2e- 99 



X 3 (0)' 



y 



X 2 (0) 



(3.50) 



(3.51) 



follows. Here, we use dimensionless variables as e = 
e/Aoo, x = x/^o and y = y/£o and we regard X{0) as 
A(fc). From eqs. (3.15)-(3.17), we obtain the parametric 
equations: 



X 
Y 



X 2 {0) 



24^cos<9-sin<9 
2^p s in(9 + cos(9 



(3.52) 



for the enveloping curve, on which the LDOS diverges. 

We consider two examples. First, we consider isotropic 
s-wave superconductor with circular Fermi surface, 
where X(0) = 1 .In this case, two paths specified by 
0\ and 02 contribute to v(r,e) in eq. (3.47), as shown 
in Fig. 2. The parametric equations for the enveloping 
curve is given by 



X 
Y 



— S\Ii0 
cos 



(3.53) 



which forms the circle with a radius of e as shown in 
Fig. 3. 

Second, we consider the d-wave superconductor with 
the circular Fermi surface; X(0) = cos2#. In this system, 




Fig. 3. Set of quasiparticle paths for e/Aoo = 0.01 in the isotropic 
s-wave superconductor. The LDOS diverges on the enveloping 
curve of the quasiparticle paths, which in the present case forms 
the circle X 2 + Y 2 = I 2 . 



we obtain the following parametric equations: 



X 
Y 



cos 2 2(9 



-4cos#tan2# - sin (9 
-4 sin # tan 2(9 + cos# 



(3.54) 



where the LDOS diverges. In Fig. 4, we show the set of 
quasiparticle paths (Fig. 4(a)) and the enveloping curve 
of them (Fig. 4(b)). We can see that the enveloping 




Fig. 4. Distribution of the LDOS in the d-wave superconductor 
about e/Aoo = 0.01. (a): the set of quasiparticle path, (b): the 
enveloping curve of the quasiparticle paths. The LDOS diverges 
on the enveloping curve of the quasiparticle paths. 



curve describes the overall features of the spatial distri- 
bution of the LDOS obtained by numerical calculations 
by Schopohl et a/. 5 and Ichioka et al. 6 The LDOS pat- 
terns for different energies reduce to a single figure after 
scaling a spatial coordinate. Therefore, the LDOS pat- 
tern in the limit e — > is expressed as 



±X. 



(3.55) 



This is the reason why the LDOS pattern rotates by 45 
degree with increasing the energy. 6 

Now we return to general case and we show that the 
LDOS has the square-root singularity near the envelop- 
ing curves. We consider the LDOS on the line running 
in the radial direction as shown Fig. 5(a). We expand 
eq. (3.48) around a divergence point A(r^,6^) at which 



the line passes the enveloping curve: 
dh 



h(r,0) 



dr 



ld 2 h 
Ar H 

(rM 2d0 2 



(AO) 2 



d 2 h 



dOdr 



a a. 19 h 
ArA(9 + -— r 

(rM 2 dr 2 



(r<A) 



(Ar) 2 



(3.56) 

Since E(r, 0) ex y = rsm(6 r — 0) has the form of rx (a 
function of 0), it follows that 

d 2 h _ 1 dh 

(ri,9i) n 86 



86dr 

from eq. (3.49) and 



(ri,0») 



0, 



(3.57) 



d 2 h 
dr 2 



0. 



Hence we obtain 
dh 
dr 



h(r,0) 



- — 



(AO)' 



(3.58) 
(3.59) 



The solution of eq. (3.48) with eq. (3.59) 



(A<9)=± -2 



dh 
dr 



d 2 h 

oe 2 



(ri,0i) 



Ar 



1/2 



exists only when 



/a ^9h 



d 2 h 



<o, 



(3.60) 



(3.61) 



is satisfied. This results from the fact that quasiparticle 
paths exist only in one side of the enveloping curve (See 
Fig. 5(a)). We similarly expand eq. (3.49) around (r^,^) 
as 



dh d 2 h 
06 ~ W 



AO H 

(nji) dOdr 



d 2 h 
Ar = — - 

(r*A) 00 2 



{r u 6i) 



A6. 



(3.62) 

From eqs.(3.47), (3.59) and (3.62), the LDOS is written 

as 



v(t, e) ex 



d 2 h 
dO 2 



dh 



(ri,6i) 



Ar 



-1/2 



(3.63) 



(nA) dr 

Therefore, the LDOS has the square-root divergence as 
shown in Fig. 5(b). 
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3.1.3 Anisotropic Fermi surface 

We consider superconductors with the anisotropic 
Fermi surface. Replacing with |uf(0)| = VFg(0) and 
using eq. (2.7), we can extend the present discussion to 
the anisotropic Fermi surface because the Eilenberger 
equation is the differential equation on the line in the 
direction of the group velocity. 23 By replacing h(0) with 
e — E(0)/g(9), we see that the points where the LDOS 
diverges are given by 



x — 



y 



dX 



2! 

X 2 ' 



(3.64) 
(3.65) 



3.2 Three-dimensional isotropic Fermi surface 

We extend our formulation to superconductors with 
the three-dimensional Fermi surface. We consider a sin- 
gle vortex along the Z axis which tilts from the crystal 
c axis by angle as shown in Fig. 6(a). When we take 



0) J" 1 
vortex 

id 


c 


















a 






Fig. 6. Frame of the space, (a): the relation between the coordi- 
nates fixed to the crystal axis (a, b and c) and the other coordi- 
nates fixed to the magnetic fields (X , Y and Z). (b): the planes 
with different x- 



the X axis on the a-b plane, the Y axis is automatically 
determined. We denote by &m, &m and cm the unit vec- 
tors along X, Y, Z axis, respectively. The pair-potential 
A(r, fc), At(r, fc) in the Riccati equation (3.8) and hence 
a(r, fc) and 6(r, fc) depend on r only through X and Y 
with 



r = Xclm + Yby[ + Zcm, 



(3.66) 



because of a translational symmetry along the Z axis (i.e. 
the direction of vortex). As a result, v • V in the Riccati 
equation can be replaced by 

d d 
v- V = v x — + v Y — , (3.67) 



with v • clm — vx and v • 6m = vy- Consequently, the 
Riccati equation for three-dimensional system turns into 
the same form as that for two-dimensional system. The 
momentum and the Fermi velocity at each point on 
the Fermi surface are written as 

fc = sin x cos Qclm + sin x sin Qbyi 

+ cosxcm, (3.68) 

= kxCLM + kybM + kzCM, (3.69) 



v-pv = vf sin x( cos O&m + sin #&m) 



• cos xcm 



(3.70) 



As a real-space coordinate on om~^m plane, we introduce 

x, y by 



r = xv + VU + Zc M , 



with 



u 



(3.71) 
(3.72) 



■ sin Oclm + cos Obyi 

The resultant matrix Riccati equations for three 
dimensional system are then written as 

da 



vpsiny— — h 2u n d + aA^a — A = 
ox 

db 

UFsinx — 2uj n b — bAb + A"*" = 

ox 



with the pair-potential 
A(r,fc) = 
A f (r,fc) = 

A(fc) = 



f(r) exp(i# r )A(fc), 
/(r)exp(-W r )At(fe), 

trA(fc)At(fc). 



0, (3.73) 
0, (3.74) 

(3.75) 
(3.76) 

(3.77) 



Here r = ^ x 2 + y 2 . For the moment, we denote A(fc) by 
A for simplicity. 

The fc-dependence of A(fc) on the Fermi surface re- 
spects the direction of crystal axes (a, 6, c). Thus, it is 
convenient to write A(fc) = A(fc a , fc&, k c ) in terms of com- 
ponents of fc for each crystal axis. Since the direction of 
magnetic field (or vortex) is tilted by (j) from the c axis 
within b-c plane, fc a , k c are, respectively, written as 



k a 


= kx = sin x cos 


(3.78) 


h 


= ky cos cj) — kz sin 0, 


(3.79) 




= sin x sin cos <p — cos x sin (/>, 


(3.80) 


kc 


= fey sin + fez cos 0, 


(3.81) 




= sin x sin sin <\> + cos x cos 0. 


(3.82) 



Therefore, the pair-potential is rewritten as 
A(fc) = A ^sin x cos 0, sin x sin cos — cos x sin 0, 



sin x sin 6/ sm < 



cos x cos 



6) .(3.83) 



Using eq. (3.83), we can obtain the LDOS as a function 
of r, e and (j). In what follows, ^-dependence of A(fc) will 
not be shown explicitly. 

Near the vortex core (|r| <C £), the local density of 
states is given by 



z/(r, e) = lim 



z/(0> F 



5^0 4^o 
^(0)^F 



dX sin x 

^0 



dO XSe~ 2XF ^ 

2k (e-E 3D ) 2 + 5 2 ' 
(3.84) 



/»7T /»Z7T 

/ dxsinx / dOX5(e - E 3 r>)t 
Jo Jo 



-2XF(x) 



(3.85) 
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\Wx)\ 



smx, 



(3.86) 



y 



-AooA 2 (e,x). (3.87) 



with the approximations 

c(e, x ) ~ 

Co smx 

Now we discuss the condition that v(r,e) diverges at 
a given point ro for an energy e. Let (0o,Xo, r o) be a 
solution satisfying 

h(0, X ,r)=e-E SD (0,x,r)=0. (3.88) 

To keep satisfying eq. (3.88), A9 = 9 — 6o, A% = x~ Xo 
and Ar = r — r have to satisfy 

dh 
89 



o dx 



dh 

o or 



■Ar = 0, 



(3.89) 



where || 
that 



9/i 



(eo,xo,r ) 



do 



and so on. First we assume 



dh 
dx 



> 



and we introduce a set of new variables 



dh 



o dx 



Ax, C="^- 
o dx 



o aw 



(3.90) 



A X - 
(3.91) 



We then obtain 



d0dx5(h(0,x,r)) = -) x2 - \ 2 . (3.92) 



f dh \ i ( 'dh \ 



The denominator in the right hand side is nonzero under 
the condition eq. (3.90), and therefore v(r,e) is finite. 
When 



dh 
06 



_ dh 

o dx 







(3.93) 



hold, Af9, Ax and Ar have to satisfy 

dh 



a(A0y + 6(Ax) 2 + cA(9Ax + — • Ar = 0, (3.94) 



<9r 



with 



_ ld 2 h 



h-- — 

o' " 2cV 



d 2 /i 



<9<9<9x 



(3.95) 



in order to satisfy eq. (3.88). In terms of a set of new 
variables, 



AO 

A X 

tan 2(p 
we obtain 
z>(r,e) 



cos (/? — sin 
sin (/? cos if 



P 
Q 



a — b J 

J d6d X S(h(6,x,r)) 

5(aP 2 + f3Q 2 + ^ 



with 

a = i(a + 6+V(«-^) 2 + c 2 ), (3.100) 

(3 = \(a + b- V(a-^) 2 + c 2 ). (3.101) 

z/(r, e) and z/(r, e) have the same singularity around 
(#0, Xo, ^o)- Therefore we consider z>(r,e), instead of 
z/(r, e). First we show that z/(r = ro, e) does not diverge 
when af3 > 0; when a > and /? > 0, £(ro, e) is written 
as 

pR c p27t 

z>(r ,e) oc 5(R 2 )RdRd$, (3.102) 

Jo Jo 

R 2 

= 7T f ° 5(R 2 )d(R 2 ) = J, (3.103) 
Jo 2 

with ^faP = Rcos^ and = i?sin<£. When aft < 0, 
on the other hand, v(r = ro, e) diverges logarithmically. 
When a > and /? < 0, z/(r, e) becomes 

^ c r2n RdRd<$> 



z>(r,e) 



// 

Jo Jo 



/~aj3 



xS{R 2 cos 2$ + ^ 
ar 



•Ar 



1 



2 A /a6- 



Jo 




d$ 



(3.104) 



x6 i? 2 cos2$- 



cos2<£ 
9ft 



<9r 



Ar 



In 



oc 



dh 
dr 



•Ar 



,(3.105) 



(3.106) 



where 6 is Heaviside step function. Note that the condi- 
tion af3 = ab — c 2 /4 < is equivalent to the saddle point 
condition 

d 2 hd 2 h ( d 2 h\ 2 n /oi _, 

Wd7 2 -{do^) <0 ' (3 ' 107) 

for a surface given by h(6,x) = const. The logarith- 
mic singularity of the LDOS near the enveloping curve 
is schematically shown in Fig. 7 



- / s ( 

JV3(0,0) \ 



dh 



dr 



(3.96) 
(3.97) 

(3.98) 
•Ar ) dPdQ, 





(3.99) 



Fig. 7. (a) Path in the radial direction. The set of the quasiclas- 
sical paths for a fixed % are drawn, (b) Local density of states on 
the path in the radial direction. The singularity is logarithmic. 



We note that the singularity of the LDOS in 
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three-dimensional system is weaker than that in two- 
dimensional system. Further, in three-dimensional sys- 
tem, the singularity of the LDOS exists in both sides of 
the enveloping curves as seen in Fig. 7(b). This fact is in 
contrast to two-dimensional system. These two features 
are understood in the following way. The LDOS in three- 
dimensional system can be regarded as the sum of the 
LDOS of fictitious two-dimensional system perpendicu- 
lar to the Z axis with reduced Fermi velocity vpsinx- 
Repeating the discussion on two-dimensional system, at 
the points where 



y 



esinxJ! 
1 A3 
esinx 

A 2 



(3.108) 
(3.109) 



are satisfied, a singular contribution comes from each 
two-dimensional system. Each two-dimensional system 
has a different Fermi velocity vpsinx- Hence, singu- 
lar contributions are gradually shifted in real space 
(Fig. 8(a)). Owing to the integration with respect to x? 
the resultant LDOS becomes less singular but the singu- 
lar behavior exists in both sides of an enveloping curve 
(Fig. 8(b)). 





Fig. 8. (a): Superposition of the divergence point of the LDOS 
with various x- The LDOS in the three-dimensional system. 



For example, we consider an isotropic s-wave super- 
conductor with spherical Fermi surface (A(0, x) = !)• I n 
this system, we obtain the parametric equations for a 
fixed x : 



X 
Y 



esmx 



— sin# 
cos# 



(3.110) 



where the LDOS would diverge if the system were two- 
dimensional. In Fig. 9, we show the superposition of the 
patterns with various x (Fig. 6(b)). From the LDOS 
patterns with various x planes, we can infer the points 
where the LDOS diverges in the three-dimensional sys- 
tem. From Fig. 9, we see that these points are on the 
circle with a radius of e in this system. 

In the case of the two-dimensional LDOS pattern on 
the plane fixed to x? we solve 



h(0,x) 



y 



smx 



■\ 2 (e,x)=o, 



(3.111) 




-0.02 -0.015 -0 



II 0.015 0.02 



Fig. 9. Superposition of the LDOS patterns with various %• 



The equations dh/dx = and eq. (3.107) can be re- 
garded as a consequence of three-dimensionality of the 
Fermi surface. These equations relate and x- For ex- 
ample, in the system where the pair-potential does not 
depend on x, dh/dx = is written as 



dh 
dx 



cosx 

e 

smx 



0, 



(3.113) 



with y = esinx/A 2 . Therefore, we can understand that 
X = 7r/2 satisfies this relation. The saddle point condition 
(3.107) is written as 



d 2 h d 2 h 



d9 2 d X 2 \d6dx) 



(3.114) 



with x = 7r/2 and y — esinx/A 2 , which is always satis- 
fied. Therefore, if the pair-potential does not depend on 
X, the points where the LDOS with three-dimensional 
Fermi surface diverges are equal to the points where the 
two-dimensional LDOS on the x = plane diverges. In 
the isotropic s- wave superconductor with spherical Fermi 
surface, the LDOS pattern is the circle with radius of e 
from eq. (3.110). If the pair-potential depends on x, the 
LDOS patten is more complicated. 

When we consider the anisotropic Fermi surface in the 
three-dimensional system, we just have to replace with 
|vf(0)| = VYg(0,x)' As in the case of two-dimensional 
system, we can apply our formulation to the case of 
anisotropic Fermi surface. 

4. Local Density of States 

In ref. 15, we showed the LDOS patterns on the 
heavy Fermion superconductor without inversion sym- 
metry CePtaSi, the calculation of which is outlined in 
Appendix B. In the present section, we show the LDOS 
patterns on the various materials such as NbSe2 and 
YNi 2 B 2 C. We use eq. (3.42) or eq. (3.84) with finite S 
to calculate the distribution of the LDOS with the nu- 
merical integration of eq. (3.43) or eq. (3.85). 



dh 
dO 



smx 



-A 2 (#,x) 



sinx 



A(0,x)=O. (3.112) 



J. Phys. Soc. Jpn. 



Full Paper 



Author Name 9 




4-1 Case of NbSe2-' singlet Cooper pairing and two- 
dimensional Fermi surface 
We assume that NbSe2 has an anisotropic s-wave pair- 
ing, which is written as 

A(r,fe) = ia y i/j(r,k), (4.115) 

^(r,fc) = |A(r)|exp(i0 r )(l + c A cos60). (4.116) 

Here, r is the phase in real space and ca denotes 
the degree of anisotropy in the superconducting energy 
gap. The experimental data of scanning tunneling spec- 
troscopy indicate that ca ~ 1/3. 4 In the magnetic field 
parallel to the c axis, we obtain A written as 



1 + \ cos 6(9. 



(4.117) 



We assume the two-dimensional Fermi surface because 
NbSe2 is a layered superconductor. We take g = 1 for 
simplicity. Figure 10(a) shows the set of quasiclassical 
paths of Andreev bound states with energy e = O.OlAoo. 
The parametric equations for the enveloping curves are 
given by 



X 



Y 



(1 + I cos 6(9) 2 



cos sin 60 
: (1 + I cos 6(9) 

, sin sin 60 



sin (9 J , (4.118) 
+ cos0 J .(4.119) 



(1 + § cos 6(9) 2 V (1 + 1 cos 6(9) 

We show the points where the LDOS diverges in 
Fig. 10(b). In Fig. 11, we show the distribution of the 
LDOS with the numerical integration of eq. (3.42). Our 




Fig. 10. Distribution of the LDOS in NbSe 2 about e/Aoo = 0.01. 
(a): the set of quasiparticle path, (b): the enveloping curve of the 
quasiparticle paths, the LDOS diverges on the enveloping curve 
of the quasiparticle paths. 



LDOS pattern is consistent with earlier numerical results 
on the LDOS pattern calculated by Hayashi et al. 4 

4.2 Case of YN12B2C: singlet Cooper pairing and 
three-dimensional Fermi surface 
Next, we consider YNi 2 B 2 C. The LDOS was probed 
by STM in this material. 2 We assume that YM2B2C has 
point nodes in quasiparticle gaps. 24 26 The pair-potential 
is written as 



A(r,fe) = ia y i()(r,k), 



(4.120) 



|A(r)|exp(*0 r )i(l 



sin 4 X cos 4(9), (4.121) 



Fig. 11. Distribtion of the local density of states about NbSe2- 
Here, e/Aoo = 0.01 and the smearing factor = 0.001. The 

horizontal axis corresponds to the a axis. 



-(1 — sin 4 x cos 40). 



(4.122) 



Here, we assume the spherical Fermi surface (g = 1) 
and we consider the magnetic field parallel to the c axis. 
From eq. (3.93), we obtain two parametric equations for 
enveloping curves because the number of the solutions of 
dh/dx — is two. One solution is given by cosx = so 
that the parametric equation is written as 



X 



4e 



(1 -cos40) 2 
4e 



(1 -cos4(9) 2 V (1 - cos 4(9) 



sin 40 cos 
'(1- cos 4(9) 

sin 40 sin 



sm( 



+ cos 



(4.123) 



(4.124) 



The other solution is given by sin 4 x = — (7 cos 4$) 1 so 
that the other parametric equation is written as 



X 



Y 



with 



49 
~16 e 

49. 
"IT 



7 cos 4(9 



7 cos 40 



< 



1/4 



(tan 4(9 cos (9 + sin 0) ,(4.125) 



1/4 



(tan 4(9 sin - cos 0) ,(4.126) 



1 



< 1. 



(4.127) 



7 cos 40 

The points where the LDOS diverges are two curves 
shown in Fig. 12. One curve runs away to infinity, while 
the other curve stays within a finite distance from the 
origin. We also calculate the distribution of the LDOS 
as shown in Fig. 13 with eq. (3.84). The only lines ap- 
proaching to infinity are experimentally observed at zero- 
energy in the case of the superconductors with gap nodes. 
In the case of YNi2B2C, for example, the cross-shaped 
LDOS pattern can be observed at zero-energy, where the 
nodes exist at (0,x) = (0,7r/2), (7r/2,7r/2), (7r,7r/2) and 
(37r/2, 7r/2) (See, Fig. 13(d)). With energy increasing, the 
cross-shaped LDOS pattern rotates by 45 degree (See, 
Fig. 13(b) and (d)). 

We find that the contribution of the equator of the 
Fermi surface (x = tt/2) to the LDOS pattern is ex- 
tremely large in superconductors with three-dimensional 
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states diverges. The pair-potential is taken to be eqs. (4.120)- to the a and c axes, respectively. 

(4.122), which is related to YNi 2 B 2 C. Here, e/Aoo = 0.01. 



(a) 



(b) 




(c) 



(d) 



* 



Fig. 13. Distribtion of the local density of states for the vortex 
parallel to the c axis about YNi2E>2C. (a) e/Aoo = 1.0 xlO -2 and 
the smearing factor 5/Aoo = 0.25xlCT 2 : (b) e/Aoo = l.OxlCT 3 
and 5/Aoo = 0.5 x lO" 3 : (c) e/Aoo = 1.0 x 10~ 4 and 5/Aoo = 
0.95 x 10- 4 : (d) e/Aoo = 1.0 x lO" 5 and 5/Aoo = 0.95 x 10~ 5 . 
The horizontal and vertical axes correspond to the a and b axes, 
respectively. 



Fermi surface. Therefore the LDOS pattern of YM2B2C 
is similar to that of two-dimensional d-wave supercon- 
ductor. By the rotation of the magnetic fields, one can 
obtain the three-dimensional infomation of the pairing 
symmetry. For example, in the magnetic field parallel to 
the a axis, the LDOS pattern of YM2B2C reflects the 
two-nodes property of the equator on the Fermi surface 
as shown in Fig. 14. If YM2B2C is a d-wave supercon- 
ductor, a line node is located along the equator and the 
LDOS looks different from that shown in Fig. 14. 



5. Remarks on the Formulation 

We expanded the coefficients of Riccati equations 
a(x, y, e) written as 



da(x, y = 0, e 



0)^ + da(x,y = 0, e = 0)^ 



dy :> de 

(5.128) 

with respect to y and e. This expansion is valid when 
\y\/£o <C 1 and |e|/Aoo <C 1 are satisfied, since typical 
length and energy are, respectively, £0 and A^. Although 
it seems that the expansion of a pair-potential written 
as eqs. (3.24) and (3.26) with respect to y is dangerous 
near x = 0, we note that the dimensionless expansion 
parameter is not y/x but y/£o- It should be noted that 
the appropriate region of our approximation is |#|/£o < 1 
because the perturbations of a and b can be applied in the 
regions — 00 < x < £0 and — £0 < x < 00, respectively. 
Practically, this condition does not matter so much since 
the LDOS decreaces exponentially as exp(— |r|/£o) away 
from the vortex center. 

In our analytical theory, there are two advantages in 
comparison with direct numerical calculations of Riccati 
equations. First, we can obtain physical interpretations 
of the LDOS pattern. The LDOS pattern consists of the 
enveloping curve of the quasiclassical paths. We can see 
from which part of the Fermi surface large LDOS at a 
point comes from. Second, we can calculate the LDOS 
in any anisotropic superconductors and in any directions 
of magnetic fields, analytically more easily than numeri- 
cally. 

We calculated the local density of states around a 
vortex core in the cases of NbSe2 and YNi2B2C. From 
eqs. (3.108) and (3.109), we see that the LDOS patterns 
are in proportion to e. The LDOS diverges at the origin 
at zero-energy. If the energy e becomes large, the size of 
the LDOS patterns becomes large. Therefore, one can ob- 
serve more detailed structure of the LDOS pattern with 
higher energy region. 

Besides the LDOS, we can also calculate the physical 
quantities expressed as an one-particle operator because 
the quasiclassical Green function around a vortex core 
is derived. For example, the distribution of the current 
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density around a vortex core is written as 

j(r) = 2ei7rT^2 J ^^Trg(k,r;iu n ). (5.129) 

n 

In real materials, the LDOS does not diverge because 
of the smearing effect. Both effects of phonons and im- 
purities contribute to the smearing effect. Therefore, the 
smearing factor S in eq. (3.42) or (3.84) determines the 
cleanness of the LDOS pattern. The smearing factor can 
be estimated from the coherence length and the mean 
free path. Since the smearing factor enters into Green 
function only through e + i<5/Aoo, the LDOS depends 
mainly on e or ^/Aoo, whichever is larger. In the low- 
energy limit e <C ^/Aoo, the LDOS patterns do not de- 
pend on 6 because of the smearing effect; zero bias ex- 
perimental data should be compared with our result in 
clean limit for energy comparable to electron-phonon re- 
laxation rate or impurity scattering rate. 

In Appendix B, we show how to solve the Riccati equa- 
tion in the case of CePtaSi in which the Fermi surface is 
split by spin-degeneracy lifting due to a spin-orbit cou- 
pling. This method has been used in ref. 15. However, 
diagonalizing the Gor'kov equation in the spin space be- 
fore the application of the quasiclassical theory, we ob- 
tain the two diagonalized Riccati equations correspond- 
ing to the split Fermi surface by spin-orbit coupling. We 
do not assume a <C A in the method diagonalizing the 
Gor'kov equation (a represents the strength of the spin- 
orbit coupling) and the resultant Green function in this 
diagonalizing method is equivalent to the Green func- 
tion obtained in Appendix B. Therefore, our theoretical 
formulation has a wider applicability. 

6. Comparison with the STM Experiment 

We compare our results with the experiment with STM 
by Nishimori et al. 2 They observed the LDOS patterns by 
STM in the isolated vortex of YM2B2C and found that 
the LDOS around a vortex core is fourfold-symmetric 
star-shaped in real space (See, e.g., Fig. 2 in ref. 2 ). 
The band width W obtained by the numerical calcula- 
tion of the band structure in YM2B2C 27 is W ~ leV and 
the gap energy is A ~ lmeV so that A/W ~ 1/1000. 
This suggests that the quasiclassical theory is appliable 
to YM2B2C. The existence of gap nodes makes the en- 
ergy levels continuous, according to a quasiclassical cal- 
culation supplemented by the Bohr- Sommerfeld- Wilson 
quantization condition. 16 Further, the result about the 
lowest energy levels of the bound state calculated by 
Hayashi et al. 28 cannot be applied to YM2B2C because 
their calculation is in the case of isotropic s-wave super- 
conductor. Therefore, it is not clear whether the STM 
experiment on YM2B2C is the observation of quantum 
regime vortex core. The size of the core structure in this 
STM experiment is comparable to the coherence length 
£o- On the other hand, the appropriate region of our ap- 
proximation is r <C £o- It might be beyond the scope 
of our theory to compare our results with experimental 
results. However, the LDOS pattern by our calculation 
about NbSe2 is consistent with the pattern of the numer- 
ical calculation by Hayashi et a/., 4 which is consistent 
with the STM experiment by Hess et al. 1 Therefore, we 



believe that it is physically meaningful to compare our 
results with experimental results. 

From Fig. 3(b) in ref. 2, it is noticed that the angu- 
lar dependence of N s (0 meV,r,#) does not depend on 
the radial distance from the core. This means that the 
arm of the star-shape pattern spreads. From eqs. (4.123)- 
(4.126), the arm of the star-shape pattern does not 
spread in the zero-energy limit. However, considering the 
smearing effect, the LDOS pattern in the zero-energy 
limit with the smearing factor S is similar to the pattern 
with the energy e = S. (Refer to the discussion in §5.) 
Therefore, we compare the LDOS shown in Fig. 13(d) 
with the experimental result. We can find that the arm of 
the LDOS pattern shown in Fig. 13(d) spreads. As shown 
in Fig. 15, the angular dependence of the LDOS does not 
depend on the radial distance from the core. Thus, the 
spreading star-shape pattern can be understood by our 
theory. In the large energy region, on the other hand, 
the location of each peak of the LDOS is different from 
that in the STM experiment (See, Fig. 16 in this paper 
and Fig. 5(d) in ref. 2). This suggests that more realistic 
approximation is necessary for more qualitative compar- 
ison with the STM experiment. For example, we have to 
calculate the LDOS with the anisotropic Fermi surface, 
for example, tight-binding model. 

0.03 1 1 1 1 1 1 1 1 1 

0.07 
0.06 - 
0.05 - / 
^0.04 - i 
0,03 i 




0.01 0,02 0.03 0.04 0.05 0,06 0.07 0.08 



Fig. 15. Radial dependence of the width of the arm of the star- 
shape LDOS pattern about YNi 2 B 2 C. e/Aoo = 1.0 x 1CT 5 and 
(5/Aqo = 0.95 x 10 -5 . The angular dependence of the LDOS does 
not depend on the radial distance from the core. 




Fig. 16. LDOS pattern in the large energy region. The intensity 
of the LDOS is large near the four circles. 

Our results show that the pairing symmetry itself can 
be investigated by the detailed observation near a sin- 
gle vortex core (r <C £o)- The past experiments by STM 
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have been made in the region r ~ £. In this region, it is 
hard to investigate rich structure of the LDOS pattern 
because of the exponentially reduction of the intensity 
away from the core. Therefore, If the LDOS in much 
narrower energy-range near zero energy and with much 
higher spatial resolution are observed, more detailed in- 
formation of the pairing symmetry will be obtained. 

7. Conclusions 

We investigated the local density of states (LDOS) 
around a vortex core on the basis of the quasiclassical 
theory of superconductivity. We obtained the analytical 
formulation of the LDOS in any anisotropic supercon- 
ductors and in any directions of magnetic fields with the 
use of Kramer-Pesch approximation. We calculated the 
LDOS pattern of NbSe 2 and YNi 2 B 2 C. The LDOS pat- 
tern of NbSe 2 is sixfold-symmetric star-shaped and con- 
sistent with the numerical calculation by Hayashi et al 4 
The LDOS pattern of YNi 2 B 2 C is fourfold-symmetric 
star-shaped near the zero-energy, which is qualitatively 
consistent with the STM experiment. Further, without 
invoking to quantum effect, we exhibited the possible 
reason why the low-energy quasiparticle states spreading 
around the a axis direction were observed in the STM ex- 
periment. The Band width W estimated from a band cal- 
culation 27 and the experimentally observed gap energy 
show that A/W - 1/1000, and therefore YNi 2 B 2 C seems 
not to be in the quantum regime. On the other hand, if 
the LDOS pattern of YNi 2 B 2 C in the magnetic field par- 
allel to the a axis is similar to Fig. 14, YNi 2 B 2 C is likely 
to be a point-node superconductor. Our results explic- 
itly demonstrate the way to determine the anisotropy of 
pairing-potential in superconductors from STM experi- 
ments. 
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Appendix A: Projection Method 

In this Appendix, we simplify the Eilenberger equation 
with a parametrization for the propagators, which were 
originally introduced in the study of the vortex dynamics 
by Eschrig. 10 ' 29 First, we define the projectors written as 

P± = \{iT'g). (A-i) 

Here, these projectors satisfy the following relations: 

P±-P± = P±, (A-2) 



P+-P- 
P++P- 



P--P+ 



6, 



(A-3) 

1, (A-4) 

g = P--P+, (A-5) 

= -(2P+ - 1) = -(1 - 2P_). (A-6) 

Substituting these projectors into eq. (2.4), we obtain the 
equation of projectors 



iv F ■ VP± + [iw„f 3 - A, P±] = 0. 



(A-7) 



Next, we introduce the matrix a and b. Using a and 6, 
we define P± as 



p+ 



i 



—a 
1 



(i-ab)-\l,a) (A-8) 



(l-Sa)- 1 ^,!), (A-9) 



which satisfies eqs. (A-2) and (A-3). To satisfy eq. (A-4), 
we define the relation between a and b, which is written 



as 



(1 — ab) 1 a = a(l — ba) 1 . 



(A-10) 



Substituting eqs. (A-8) and (A-9) into eq. (A-7), we ob- 
tain the matrix Riccati equations. 

Appendix B: Case of CePt 3 Si: Mixed Singlet- 
Triplet Cooper Pairing and Three- 
dimensional Fermi Surface 

We calculated the LDOS in the case of CePtaSi in 
ref. 15. CePtaSi is the heavy fermion superconductor 
without inversion symmetry, 30 ' 31 where the singlet pair- 
ing and the triplet pairing may be mixed. 32 ' 33 The spin- 
orbit coupling term exists in the Hamiltonian. We need 
to extend our formulation in the body of this paper be- 
cause the pairing symmetry in this material does not 
satisfy the relation AA^ oc &q. The LDOS patterns in 
CePtaSi are derived from the following discussion. This 
discussion can also be applied to the non-unitary triplet 
superconductors. 

Considering the result of numerical calculations by 
Hayashi et al, M we assume that the spatial varia- 
tions of s-wave pairing component of pair-potential and 
that of p-wave pairing component are the same: A = 



ky& X 



knr, & 11 



[Va + d k • &] ia y = A(r) 
with the s- wave pairing component Here, the ratio 
of the singlet to triplet component is defined in the bulk 
region as \j/ = ^ / A. This mixed s + p-wave model is 
proposed for CePtaSi. 33 The Eilenberger equation which 
includes the spin-orbit coupling term is given as 20 ' 35 38 

iv F (fe) • V£ + [iu; n f 3 - ag k • S - &,g] = 0, (B-l) 



with 

9k = 







9k 





9-k°o 



dk^o 






-9k°o 



(B-2) 



We neglect the impurity effect and the vector potential 
because CePtaSi is a clean and extreme- type-II super- 
conductor. 30 The Riccati equations are written as 

vy • Va+ + 2u n aj r + a+A^a^ — A 

+i(a+a(^-(T) tr + a^-(ra + ) = 0, (B-4) 

v F • VL - 2u n b- - b- AL + A f 

-i{b-ag k .& + a{g k .&)HJ) = 0. (B-5) 

With the Kramer-Pesch approximation as with the 
means in the previous sections, we expand Riccati equa- 
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tions by y, uj n and a. In the 0-th order, equations (B-4) 
and (B-5) are the same as ones in the body of this pa- 
per. To obtain the solution of these equations, we assume 
that the 0-th solutions do and bo are written as 



ao 



an 
Q>21 



a>i2 

«22 



, ^0 



&21 ^22 



(B-6) 



and substitute these solutions into the 0-th order equa- 
tions (3.28) and (3.29) with v F — > vpsinx- Solving these 
equations as simultaneous equations about four vari- 
ables, we obtain the solutions written as 



d = ± 




1 


1 



-1 















, (B-7) 



(B-8) 



with k± = k x ± iky. Here, we determine the signs of 
ao and bo by the fact that the 0-th order Green function 
go = — (l+ao6o) _1 (l— dobo) diverges (See §3.1.2), namely 
&o^o = — 1- We consider the limit of |rf| — > and that 
of \£ — > 0, where pair-potentials are unitary. In the limit 
of \d\ — > 0, the pair-potential has the singlet component 
only. From eq. (3.30), ao is equal to ia y . In the limit of 
\£ — > 0, the pair-potential has only the triplet component. 
In the first order, eqs. (B-4) and (B-5) are written as 



v F smx 



ddi 
dx 



v F sin x 



dx 



-f 2u n d + a AJai + aiAJa + a A{a 
Ai + i(a a(g k • <x) tr + ag k • <ja ) = 0, (B-9) 
2uj n bo - boAoh - &iA &o - 6 Ai6 



+A f x - i(b ag k - <t + a(g k - <r) tr 6 ) = 0. (B-10) 

To solve the above equations, we consider the solution of 
the homogeneous equations written as 



ddi 



v F sin + a AJai + d 1 A J d 

dbi t X t i X i 

v F smx- b A bi - DiA D 

ox 



vt; 



We assume the solution: 

ai = G exp 

V v F sin x 



/ d^A^) 



0, (B-ll) 
0. (B-12) 

(B-13) 



and substitute this into eq. (B-ll). Therefore we regard 
the equations as the equations with variables of vector 
A = (a, 6, c, d) and find the eigenvalues (= —A) and 
eigenvectors Ai. The same procedure is applied to b\. 
With the method of variation of constants [C — > C(x)], 
namely with writing 

4 / cw Ax\ \ 



di = y^Cj(x) exp 



2A ?; 



2=1 



smx Jo 



1 



dx'A(x') )Ai, (B-14) 



eqs. (B-9) and (B-10) are solved, so that we obtain the 
solutions of eqs. (B-4) and (B-5) up to the first order. 
The explicit expressions for the resultant Green function 

and the LDOS pattern are given in ref. 15. 
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